AD-A213  145 


Wr?- 


TECH.  MEMO 
SPACE  370 


UNLIMITED 


b }  \  1  I  0  >  7  >  (^) 


TECH.  MEMO 
SPACE  370 


ROYAL  AEROSPACE  ESTABLISHMENT 


I  ' ON  UNIVERSAL  ELEMENTS,  AND  CONVERSION  PROCEDURES  TO  .s;'' 

7;:^  -V  •/-/  .....  :  AND  FROM  POSITION  AND  VELOCITY  : 


■  vv>w: 

A  -li.-iTU*'  - 


gfiDTIC 

ELECTE  S 
^  OCT  0  6  1989 1 


R.  H.  Gooding 


:rj> .  :sr  ^:.Tir»i<55*2S  ^  *-»r-5r  ■£: 

I  .  ■.•’i  L'r'AY  Ai 


!  — .  -  '  ^  -  *-tj  V  •  -» 


iriSSrifias^r  « v; *•  •  •%  !/*i.  iLJr  ..-"i  .  -Lr 

';■  f  ttSte?*  -  x  V-.  r  r  _i- .  ,  M  m  ~.  i  ■  ~  v  -  ■  ...  . .  .•■*  ■ .  7.  -  -■Vjit.  -V. 

»aasr#*.  Wt'-WS1*.;" -■*  - «  •.-  -■■.  ■•  ■'  .  ...  ,  ■■  .''  *  ~  :V';y 


D 


=rT-  V^K^'ggS^U 


DISTRIBUTION  STATEMEKT  A 

Appxoved  ioi  public  lelactel  , 
Diamounon  Uniimtad 


Procurement  Executive,  Ministry  of  Defence 
Farnborough,  Hants 


UNLIMITED 


89  10  6  029 


0047263 


CONDITIONS  OF  RELEASE 


BR-  11099t 


COPYRIGHT  (c) 
<983 

CONTROLLER 
HMSO LONDON 


********* . *******  Y 

Reports  quoted  are  not  necessarily  available  to  members  ol  the  public  or  to  commercial 
organisations. 


UNLIMITED 


ROYAL  AEROSPACE  ESTABLISHMENT 

Technical  Memorandum  Space  370 
Received  for  printing  4  July  1989 


ON  UNIVERSAL  ELEMENTS,  AND  CONVERSION  PROCEDURES  TO 
AND  FROM  POSITION  AND  VELOCITY 

by 

R.  H.  Gooding 


SUMMARY 


An  element  set  is  advocated  that  is  familiar  (in  traditional  terms),  and 
yet  applicable  to  all  types  of  orbit  without  loss  of  accuracy.  It  is  not  free  of 
singularity,  but  this  is  not  a  serious  deficiency.  Conversion  procedures,  to  and 
from  position  and  velocity,  are  outlined,  with  Fortran-77  listings  appended. 

Tests  have  indicated  that  the  errors  in  the  pair  of  procedures  are  minimal, 
accuracy  being  limited  only  by  computer  precision  and  the  (fixed)  number  of 
iterations  used  in  the  Kepler-equation  solutions. 

This  is  the  original  text  of  a  paper  that  has  now  been  published  in  the 
journal  Celestial  Mechanics  (Vol.44,  pp  283-298,  1988).  The  paper  is  printed 
here,  from  page  3,  in  the  format  required  by  the  journal,  the  contents  being 
listed  on  page  2.  The  paper  is  a  shortened  version  of  RAE  Technical  Report 
87043,  and  a  companion  paper,  also  published  in  Celestial  Mechanics  (Vol.44, 
pp  267-282,  1988)  is  available  as  Technical  Memorandum  Space  369. 


Copyright 

© 

Controller  HMSO  London 


1989 

H-l 


UNLIMITED 


2 


LIST  OF  CONTENTS 


1  INTRODUCTION 

2  A  SET  OF  UNIVERSALLY  APPLICABLE  ELEMENTS 

3  CONVERSION  OF  ELEMENTS  TO  POSITION  AND  VELOCITY 

4  CONVERSION  OF  POSITION  AND  VELOCITY  TO  ELEMENTS 

5  TESTING  OF  THE  COMPUTING  PROCEDURES 

6  ON  THE  SINGULARITIES 

7  CONCLUSION 

Appendix  A  Subroutine  ELS2PV 
Appendix  B  Subroutine  ELS3PV 
Appendix  C  Subroutine  PV2ELS 
Appendix  D  Subroutine  PV3ELS 
References 
Illustration 

Report  documentation  page 


Page 

3 

5 

7 

10 

13 

16 

18 

21 

22 

23 

25 

26 

Figure  1 

inside  back  cover 


TM  Sp  370 


3 


Abstract.  An  element  set  is  advocated  that  is  familiar  (in  traditional 
terms),  and  yet  applicable  to  all  types  of  orbit  without  loss  of 
accuracy.  It  is  not  free  of  singularity,  but  this  is  not  a  serious 
deficiency.  Conversion  procedures,  to  and  from  position  and  velocity, 
are  outlined,  with  Fortran-77  listings  appended.  Tests  have  indicated 
that  the  errors  in  the  pair  of  procedures  are  minimal,  accuracy  being 
limited  only  by  computer  precision  and  the  (fixed)  number  of  iterations 
used  in  the  Kepler-equation  solutions. 

1 .  Introduction 

There  has  long  been  the  goal,  in  celestial  mechanics,  of  subsuming  the 
solutions  of  the  two-body  problem,  particular  to  the  ellipse,  parabola 
and  hyperbola,  in  a  universal  solution,  valid  for  the  rectilinear  orbit 
of  each  type  as  well  as  the  general  orbit.  Steps  towards  this  goal  have 
been  taken  by  Sundman  (1912),  Stumpff  (1947),  Goodyear  (1965),  Herrick 
(1965),  Pitkin  (1965)  and  Shepperd  (1985),  among  others,  and  some 
elegant  mathematical  formulations  have  been  achieved.  The  aesthetic 
attractions  of  the  universal  approach  have  to  be  weighed  against  a 
numoer  of  computational  disadvantages,  however,  and  the  latter  are  often 
disregarded. 

An  important  distinction  must  be  made  between  universality 
associated  with  orbital  elements,  on  the  one  hand,  and  the  formulae  and 
working  variables  of  computing  procedures,  on  the  other.  The  case  *or 
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universal  elements  is  undeniable,  but  it  does  not  automatically  extend 
to  algorithms  such  as  those  involved  in  the  conversion  of  elements  to 
and  from  position  and  velocity;  the  user  of  such  an  algorithm  is 
naturally  concerned  with  its  generality,  as  well  as  with  its  accuracy 
and  efficiency,  but  a  requirement  for  generality  does  not  imply  that  the 
algorithm  has  to  be  internally  'blind'  to  orbit  type. 

The  foregoing  distinction  is  at  the  heart  of  the  present  paper, 
which  is  a  shortened  version  of  a  recent  RAE  Report  by  Gooding  (1987). 
Section  2  discusses  what  is  actually  meant  by  a  set  of  universal  (or 
universally  applicable)  elements,  and  introduces  the  set  to  be  used  in 
the  rest  of  the  paper;  this  set  suffers  from  certain  singularities,  but 
the  singularities  are  found  to  lead  to  no  real  difficulty  in  the  con¬ 
version  algorithms.  Based  on  the  assumed  element  set.  Section  3 
describes  a  computing  procedure  for  the  conversion  to  position  and  velo¬ 
city;  this  involves,  in  particular,  separate  algorithms  for  solving 
Kepler's  equation  and  the  corresponding  hyperbolic  equation,  as 
currently  described  by  Gooding  and  Odell  (1989).  Section  4  describes 
the  reverse  procedure;  here  the  algorithms  for  ellipse  and  hyperbola 
need  to  be  formulated  with  great  care,  since  serious  inaccuracy  can 
arise  if  a  formula  for  hyperbolic  orbits  is  based  on  direct  trans¬ 
mutation  of  the  corresponding  elliptic  formula.  Section  5  of  the  paper 
is  concerned  with  the  testing  and  performance  of  the  two  conversion 
procedures,  each  of  them  having  been  implemented  by  a  pair  of  Fortran-77 
subroutines,  one  to  cover  the  two-dimensional  (in-plane)  part  of  the 
conversion  and  the  other  the  three-dimensional  aspects. 

There  are  situations  in  which  the  orbital-element  singularities 
present  difficulties  that  are  more  real  than  in  the  conversion 
algorithms.  Section  6  indicates  how  such  difficulties  may  be  dealt 


with. 


2.  A  Sec  of  Universally  Applicable  Elements 

We  seek  to  define  a  set  of  universally  applicable  elements  for  motion  in 
unperturbed  orbits  about  a  centre  of  Newtonian  attraction  of  strength 
p  .  We  restrict  our  attention  to  non-redundant  element  sets,  for 
convenience  denoting  by  z  the  sextuple  of  quantities  that  such  a  set 
must  comprise;  further,  we  make  the  usual  assumption  that  the  first  five 
elements  of  z  define  the  orbital  path,  whilst  the  sixth  is  a  phase 
specifier.  Since  our  main  concern  is  with  computing  procedures  for 
converting  from  orbital  elements  to  position  and  velocity,  and  vice 
versa,  we  let  x  denote  the  sextuple  of  components  of  position  and 
velocity  in  a  convenient  system  of  rectangular  coordinates,  with  origin 
at  the  attraction  centre;  we  introduce  f  as  the  function  converting 
from  z  to  x  ,  so  that 

x  -  fU)  .  (1) 

We  shall  also  refer  to  the  Jacobian  (partial-derivative)  matrix  of  x 
with  respect  to  Z  »  and  denote  it  by  ,2  . 

If  a  particular  element  set  can  be  chosen  that  covers  every  type 
of  orbit,  then  in  principle  we  regard  these  elements  as  universal.  It 
is  implied  that  the  function  f  is  surjective,  with  range  covering  all 
possible  x  ,  but  this  is  not  enough.  We  also  require  that  J  is 
defined  (exists)  over  the  domain  of  valid  Z  ,  with  no  occurrence  of 
discontinuity  or  infinity. 

Ideally,  f  would  be  injective  as  well  as  surjective,  so  that  it 
would  have  a  unique  inverse,  with  the  matrix  J  never  singular.  Since 
non-singularity  seems  to  be  incompatible  with  universality,  however,  we 
must  give  up  the  extra  requirement,  but  for  many  purposes  (in  solving 
Lambert's  problem,  for  example)  this  is  of  little  consequence.  From  the 
surjective  property,  we  can  always  define  a  (non-unique)  function, 
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f  ^  ,  such  that  the  composition  of  f  ^  followed  by  f  is  the  identity 
over  x-space  ,  and  this  is  enough  to  demand  of  the  function  f  1  . 
Looking  ahead  to  Section  5,  however,  we  can  express  this  relation  more 
conveniently  (and  symmetrically)  as  a  property  based  on  £-space  ,  viz 

f  f"1  f  -  f  .  (2) 

We  now  proceed  to  the  identification  of  a  particular  set  of 
universally  applicable  elements.  Our  starting  point  is  the  set 
(a,  e,  i,  fl,  u,  M)  of  familiar  elements  used  for  elliptic  orbits.  The 
elements  i  ,  n  and  hi  ,  which  operate  as  Euler  angles  relating  the 
orbital  plane  to  the  chosen  axis  system,  are  already  universal,  applying 
as  well  to  parabolas  and  hyperbolas  as  to  ellipses.  The  element  M  is 
manifestly  not  universal,  on  the  other  hand,  since  it  is  identically 
zero  in  the  parabolic  limit.  It  becomes  universal  on  division  by  n 
(mean  motion),  however,  the  result  being  a  new  element  (time  from  peri- 
focus)  that  we  denote  by  t  ;  clearly,  x  is  just  the  negative  of  the 
element  T  ,  traditionally  used  as  an  alternative  to  M  .  The  four 
elements  so  far  considered  are  all  involved  in  indetermlnacies  (non¬ 
uniqueness)  due  to  singularity  (and  the  associated  non-injectivity  of 
the  function  f  ),  the  sources  of  singularity  being  rectilinear  orbits, 
circular  orbits,  and  orbits  for  which  i  -  0  or  i  «  u  . 

It  just  remains  to  consider  the  first  two  of  the  original 

elements,  and  the  combination  of  a  and  e  cannot  be  universal,  since 

a  parabola's  perifocal  distance,  given  by  q  -  a(l  -  e)  ,  and  parameter, 

2 

given  by  p  *  a(l  -  e  )  ,  would  not  then  be  defined.  If  we  replace  e 
by  q  ,  however,  the  new  pair  of  elements  (a  and  q)  would  be  universal 
if  we  allowed  a  to  be  infinite.  To  avoid  this  difficulty,  we  replace 
a  by  a  ,  where 
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ct  “  u/a  ,  (3) 

and  then  the  set  (a,  q,  i,  fl,  u,  t)  is  a  universal  set.  It  is  noted 
that  the  definition  of  a  by  (3)  maintains  the  desirable  feature  cf  the 
original  element  set  that  makes  one  element  a  function  only  of  energy; 
indeed  a  is  twice  the  negative  energy  per  unit  mass,  since 

-  2p/r  »  -a  (4) 

The  incorporation  of  the  factor  p  in  (3)  leads  to  certain  simplifi¬ 
cations  in  computing;  as  can  be  seen  from  (4),  it  also  permits  the 
element  a  to  be  well  defined  in  the  ultimate  limiting  case  in  which 
p  *  0  .  It  would  not  be  satisfactory  to  replace  q  by  p  ,  since 
(because  3q/3p  -  l/2e  ,  for  a  fixed  value  of  a)  a  circular  orbit  would 
then  involve  infinities  in  J  . 

In  terms  of  a  and  q  ,  the  following  complete  classification  of 
orbits  becomes  possible: 

a  >  0  a  ■  0  a  <  0 

q  >  0  general  ellipse  general  parabola  general  hyperbola 

q  ~  0  rectilinear  ellipse  rectilinear  parabola  rectilinear  hyperbola 

Figure  1  shows  how  q  varies  with  a/p  (  =  1/a)  for  selected  values  of  e 
The  curves  are  all  rectangular  hyperbolas,  and  it  is  easily  seen  how 
rectilinear  orbits  (points  on  the  horizontal  axis)  may  be  distinguished 
from  parabolic  orbits  (points  on  the  vertical  axis),  though  e  *  1  in 
both  cases. 


3.  Conversion  of  Elements  Co  Position  and  Velocity 
The  procedure  to  be  described  converts  from  the  assumed  element  set 
(a,  q,  i,  fl,  a) ,  t)  to  position  and  velocity  (x,  y,  z,  x,  y,  z).  It  has 
been  implemented  by  a  pair  of  Fortran-77  subroutines,  listed  in 
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Appendices  A  and  B.  Subroutine  ELS2PV  (Appendix  A)  covers  the  purely 

two-dimensional  (iu-plane)  part  of  the  conversion,  in  which  the  polar 

coordinates  r  and  u  (u  being  ’argument  of  latitude',  assuming  the 

normal  interpretation  of  the  axis  system),  together  with  the  radial  and 

transverse  velocity  components  V  and  V  ,  are  derived  from  the 

R  T 

elements  a,  q,  u>  and  t  ,  plus  u  .  The  definition  of  VT  is  such 
that  V  >  0  ;  also  it  is  assumed  that  r  >  0  ,  avoiding  the  polar- 
coordinate  singularity  at  r  -  0  ,  though  ELS2PV  operates  accurately  for 
all  points  of  a  rectilinear  orbit  other  than  the  centre  of  attraction 
itself.  Subroutine  ELS3PV  (Appendix  B)  first  calls  ELS2PV,  and  then 
uses  the  angles  u,  i  and  0  to  perform  the  appropriate  rotations  to 
complete  the  conversion. 

Details  of  the  computing  procedure  will  be  omitted  here.  The 
relevant  formulae  have  been  given,  without  proofs,  by  Gooding  (1987), 
and  further  detail  can  be  obtained  from  text-books.  Some  comments  on 
particular  aspects  of  the  procedure  are  called  for*  however. 

First,  in  much  of  ELS2PV,  different  algorithms  are  used  for  the 
ellipse,  hyperbola  and  parabola,  the  control  parameter  being  (the  sign 
of)  a  .  The  function  procedures  EKEPL  and  SHKEPL  are  used  to  solve 
Kepler's  equation  and  a  reformulated  hyperbolic  equation,  the  latter 
as  described  by  Gooding  and  Odell  (1987);  the  guaranteed  accuracies 
(relative  truncation  error)  for  these  procedures,  in  two  Iterations  of  a 
quartic  convergence  process,  are  14  decimal  digits  for  E  (eccentric 
anomaly)  and  20  digits  for  S  (*  sinh  H  ,  where  H  is  the  hyperbolic 
anomaly).  For  the  parabola,  Barker's  equation,  in  the  form 


,3  ,  .  ,  2 

d  +  6vqd  *  6y  t  , 


(5) 


is  solved  directly  for  d  ,  where  d  *  tr  .  Care  has  been  taken  to 
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obtain  maximum  accuracy  in  the  solution  of  this  cubic  equation,  in 
particular  by  the  use  of  a  refined  cube  root. 

The  next  point  concerns  the  value  of  v  (true  anomaly),  computed 
as  a  working  variable  of  ELS2PV.  For  hyperbolas  and  parabolas,  we  must 
have  |  v  |  <  ir  ,  so  that  the  subroutine  can  accommodate  an  arbitrary 
value  of  w  unambiguously,  with  included  multiples  of  2ir  carried 
directly  into  u  .  For  ellipses,  on  the  other  hand,  an  arbitrary  value 
of  v  can  arise,  reflecting  the  number  of  orbital  periods  present  in 
t  ,  so  that  the  mapping  (u>,  t)  -*■  u  is  no  longer  one-one.  This  is  an 
aspect  of  the  non-in jective  nature  of  the  function  f  of  Section  2, 
quite  apart  from  the  singularity  issue,  and  it  will  be  referred  to 
again  in  Section  A. 

As  a  third  point,  we  note  that,  for  the  ellipse,  r  and  v  are 
computed  from  formulae  in  sin  ^E  and  cos  ^E  ,  but  that  the  corre¬ 
sponding  formulae  (in  sinh  and  cosh  ^H)  are  not  available  for 
the  hyperbola,  as  it  is  S  (  -  sinh  H)  that  is  provided  by  SHKEPL. 

There  is  no  difficulty,  however,  as  both  r  and  v  can  be  expressed  in 

terms  of  C  -  1  ,  where  C  ■  cosh  H  ;  to  minimize  rounding  error,  we 
compute 

C  -  /(I  +  S2)  , 

followed  by 

C  -  1  -  S2/(C  +  1)  . 

A  final  comment  may  be  helpful  in  interpreting  the  listings  in  the 
Appendices  (it  applies  in  Section  A  as  well  as  here).  As  a  corollary  of 

the  appearance  of  p  in  definition  (3),  an  appropriate  (integral)  power 

of  p  is  naturally  associated  with  other  quantities  also;  thus  the 
Fortran  variables  E  and  M  refer  to  pe  and  pM  ,  respectively,  not 
just  e  and  M  • 


IQ 


4.  Conversion  of  Position  and  Velocity  to  Elements 

We  are  now  concerned  with  the  procedure  for  converting  from  position  and 

velocity  (x,  y,  z,  x,  y,  z)  to  the  element  set  (a,  q,  i>  ft,  w,  t ). 

Whereas  the  procedure  described  in  Section  3  constitutes  the  unambiguous 

function  f  of  Section  2,  the  procedure  to  be  described  here  merely 

defines  a  particular  inverse  function,  f  1  ,  of  the  many  satisfying 

(2).  It  has  been  implemented  by  a  pair  of  Fortran-77  subroutines, 

listed  in  Appendices  C  and  D,  that  are  the  natural  counterparts  of  those 

from  Sect?-m3.  Subroutine  PV2ELS  (Appendix  C)  covers  the  purely  two- 

dimensional  part  of  the  conversion,  computing  a,  q,  w  and  T  from 

r,  u,  V  and  V  (together  with  p),  whilst  PV3ELS  (Appendix  D)  imple- 
R  T 

ments  the  overall  procedure,  deriving  i  and  ft  as  well  as  the  input 
required  for  PV2ELS ,  which  it  calls. 

Mainly  because  of  the  semi-arbitrary  nature  of  the  procedure,  the 
f  *  subroutines  are  a  good  deal  more  complicated  than  the  f  sub¬ 
routines,  and  more  extensive  commentary  is  necessary  -  there  is  one 
simplification,  on  the  other  hand,  since  there  Is  no  longer  a  transcen¬ 
dental  equation  to  be  solved.  The  first  comment  refers  to  the  control 
parameter,  in  PV2ELS ,  for  the  type  of  orbit.  This  is  again  a  ,  given  by 
(4).  To  avoid  arbitrariness,  the  orbit  is  only  deemed  to  be  parabolic 
if  a  is  exactly  zero,  and  this  condition  is  very  unlikely  to  occur  in 
practice;  because  of  rounding  error,  this  is  true  even  when  f 
follows  an  f  that  operated  on  an  exactly  parabolic  orbit.  There  is  no 
lessening  of  accuracy  in  the  handling  of  near-parabolic  ellipses  and 
hyperbolas,  however,  so  it  is  only  the  efficiency  associated  with  the 
parabolic  formulae  that  is  lost;  if  this  were  considered  important,  then 
PV2ELS  could  be  modified  to  decree  'parabolic  orbit'  whenever  J  ra/p  | 


is  less  than  some  suitable  criterion  value. 


A  similar  situation  (still  in  PV2ELS)  applies  to  the  recognition 
of  a  circular  orbit,  which  only  happens  if  an  exactly  zero  value  of  e 
(which  is  just  a  working  variable)  is  computed.  This  special  case  has 
to  be  covered  because  it  constitutes  a  singularity,  with  the  values  of 
u>  and  t  both  becoming  indeterminate;  PV2ELS  sets  t  and  v  (another 
working  variable)  arbitrarily  to  zero,  after  which  o>  is  set  to  u  -  v 
as  with  every  other  type  of  orbit.  There  is  a  degree  of  arbitrariness 
associated  with  all  elliptic  orbits,  of  course,  due  to  the  non-injective 
nature  of  the  mapping  (<d,t)  u  ,  as  remarked  in  Section  3.  The  two 
obvious  options,  for  PV2ELS,  were  to  select  for  minimum  ]  t  |  or  minimum 
|  a)  |  .  The  former  is  the  better  choice,  because  it  reduces  rounding 

error  in  certain  circumstances,  in  particular  for  near-rectilinear 
ellipses,  and  this  option  has  been  implemented  in  the  subroutine;  the 
other  option  can  be  obtained,  however,  by  changing  the  value  of  a 
built-in  logical  variable  (L). 

It  was  remarked  in  Section  1  that  the  use  of  directly  parallel 
formulae,  for  the  ellipse  and  hyperbola,  can  lead  to  inaccuracy.  The 
reference  was  to  rounding  error,  and  a  good  example  of  such  error  arises 
with  the  computation  of  e  ,  required  so  that  q  can  be  computed  from 

q  -  p/(l  +  e)  .  (6) 

If  we  define 

c  ”  rV2/y  -  1 

a:-.!4 

s  -  r  V  W*/p  , 

R 

tnen  the  elliptic  interpretation  is  that  c  -  e  cos  E  and 

2  2  2 

s  -  e  sin  E  ,  so  that  e  is  given  by  e  -  c  +  s  .  For  the  hyperbola, 
however,  the  interpretations  of  c  and  s  are  as  e  cosh  H  and 


2  2  2 

e  ainh  H  ,  so  that  the  parallel  formula  is  e  »  c  -  s  .  The  potential 

sensitivity  to  rounding  error,  when  |  H  |  is  large,  is  obvious,  yet  this 

is  a  formula  frequently  given  in  text-books.  The  formula  appropriate 

2 

for  hyperbolas  is  simply  e  -  1  -  ap/p  ,  where  (since  a  is  negative) 
it  is  really  an  addition,  not  a  subtraction,  that  is  implied.  This 
formula  would,  in  turn,  be  inappropriate  for  ellipses  {of  the  remark  in 
Section  2  that  a  and  p  ,  as  opposed  to  a  and  q  ,  would  be  an 
unsatisfactory  pair  of  elements  for  orbits  approaching  circularity),  and 
PV2ELS  uses  the  optimum  formula  for  each  case. 

The  remaining  comments  refer  to  the  computation  of  the  quantities 
i,  ft  and  u  by  subroutine  PV3ELS,  prior  to  the  calling  of  PV2ELS.  The 
main  problem  is  singularity  again,  though  (as  in  PV2ELS)  there  is  the 
additional  source  of  arbitrariness  associated,  in  particular,  with  the 
convention  that  the  value  determined  for  ft  should  always  satisfy 
-ir  <  ft  <  ir  . 

The  predominant  singularity  relates  to  rectilinear  orbits,  and  it 
is  deemed  to  arise  only  when  the  angular-momentum  vector  is  exactly 
zero.  As  no  'orbital  plane*  is  then  defined,  it  was  necessary,  for  the 
operation  of  PV3ELS,  to  fabricate  one,  and  it  was  decided  to  do  this,  in 
principle,  by  taking  the  plane  that  contains  the  orbit  and  for  which 
i  ”  Vr  •  For  an  'axial  orbit',  perpendicular  to  the  reference  plane, 
there  is  a  subsidiary  singularity,  however,  with  the  general  principle 
inadequate  to  define  a  unique  plane;  for  definiteness  in  this  case,  the 
value  of  ft  is  decreed  to  be  zero. 

For  the  general  (non-rectilinear)  orbit,  superficially  the  only 
difficulty  arises  with  the  singularity  that  occurs  when  the  orbital 
plane  exactly  coincides  with  the  reference  plane,  and  this  is  easily 
dealt  with,  again  by  decreeing  ft  to  be  zero  -  this  is  just  like 
setting  u  and  t  to  zero  (in  PV2ELS)  for  an  exactly  circular  orbit. 
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But  there  is  another  difficulty;  it  arises  with  near-rectilinear  orbits, 
for  which  the  rounding  error  in  computing  the  components  of  the  angular- 
momentum  vector,  h  ,  can  be  serious.  The  problem  here  disappears  at 
once,  however,  if  instead  of  just  computing  h  (  ■  r  x  r)  we  compute 

r  jj  “  £  x  [  (£  x  £)  x  5]  ; 

further  detail  has  been  given  by  Gooding  (1987). 

5.  Testing  of  the  Computing  Procedures 

The  basic  philosophy  for  testing  was  that  the  two  conversion  procedures, 
as  implemented  by  the  subroutines  ELS3PV  and  PV3ELS,  would  be  used  to 
test  each  other.  The  validity  of  this  philosophy  emanated  from  two 
considerations:  first,  that  the  subroutines  were  essentially  inde¬ 

pendent,  with  (for  example)  Kepler's  equation  only  arising  in  (a  sub¬ 
routine  subordinate  to)  ELS3PV  and  rectilinear  orbits  only  having  'to 
be  recognized  in  PV3ELS;  secondly,  that  advantage  could  be  taken  of  the 
different  formulae  used  for  different  types  of  orbit,  to  make  a  careful 
study  of  continuity  across  the  transition  lines.  The  obvious  property 
to  test,  in  the  notation  of  Section  2,  was  that  ff-*  (the  composition 
of  f  1  followed  by  f)  should  be  the  identity  over  x-space  ,  with  no 
corresponding  requirement  for  f  ^f  ;  but  the  need  to  make  the  testing 
systematic  and  efficient  pointed  to  the  use  of  input  from  £-space 
rather  than  x-space  ,  and  this  is  why  the  testing  was  actually  based  on 
the  property  specified  by  Equation  (2).  The  testing  has  been  restricted 
to  the  verification  of  (2),  to  within  tolerable  rounding  error,  at  fixed 
instants  in  time,  but  this  restriction  does  not  limit  the  efficacy  of 
the  subroutines  or  the  testing,  since  to  proceed  from  x  at  cq  *  say» 
to  x  at  tj  ,  we  merely  have  to  increment  the  T -component  of 
£q  -  f  1(xQ)  ,  by  tx  -  tQ  ,  before  deriving  x1  from  f(£L)* 


To  attach  a  meaning  to  'tolerable  rounding  error'  ,  we  separate  the 

sextuple  x  into  the  pair  of  vectors  r  and  r  .  Also,  for  each 

'input'  x  (derived  from  an  actual  input  C  ),  we  write  ff  *(;«)  as 

x  +  Ax  .  Then  the  relative  errors  in  the  final  position  and  velocity 

are  taken  to  be  |  Ar  |  /r  and  |  Ar  |  /V  ,  respectively,  where  r  «  |  r  | 

and  V  »  |  i  |  •  It  might  be  hoped  that  the  ff-*  operation  would  not 

produce  relative  errors  more  than  about  a  (decimal)  order  of  magnitude 

greater  than  the  limiting  precision  of  the  computer  used;  since  this  was 
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a  PRIME  750,  of  limiting  accuracy  (double  precision)  about  10  ,  it 

was  therefore  reasonable  to  look  for  a  maximum  relative  error  no  greater 
than  2  x  10  .  This  goal  has  been  met  for  |  Ar  |  /r  without 

difficulty,  but  for  |  At  |  /V  only  in  a  modified  manner.  Before 
explaining  this,  we  remark  on  the  range  of  test  cases  covered. 

There  are  seven  degrees  of  freedom  associated  with  the  input  to 
ELS3PV,  since  y  is  an  argument  as  well  as  the  sextuple  t  .  It  was 
obviously  impracticable  to  carry  out  exhaustive  tests  over  a  seven¬ 
dimensional  space,  so  most  of  the  testing  concentrated  on  the  in-plane 
algorithms,  using  fixed  values  of  i,  and  u  .  In  these  tests,  the 
quantities  y,  a,  q  and  t  had  to  be  covered,  but  it  was  not 
necessary  to  vary  them  all  independently,  if  the  assumption  was  made 
that  orbits  of  the  same  shape  produce,  at  corresponding  points,  relative 
errors  of  the  same  order  of  magnitude.  Thus,  orbits  having  the  same 

value  of  aq/y  may  be  regarded  as  equivalent  in  shape,  this  being  the 

3  2  2 

value  of  1  -  e  ,  whilst  points  with  the  same  value  of  a  x  /y  may  be 

2 

regarded  as  in  correspondence,  this  being  the  value  of  M  for  elliptic 
orbits.  On  this  basis,  it  was  legitimate,  for  non-rectllinear  orbits,  to 
test  with  fixed  values  of  y  and  q  ,  taken  (arbitrarily)  as  64  and  1 
respectively;  these  tests,  with  two  degrees  of  freedom  (a  and  x), 
then  had  to  be  supplemented  by  one-degree  tests  (on  a  only)  for 


15 


rectilinear  orbits,  with  q  now  set  to  zero  and  t  given  a  fixed  value. 

A  wide  range  of  values  for  a  and  t  was  covered:  a  varied  from  an 

extreme  negative  value  of  -10  to  its  maximum  permissible  value  of  64 

(for  the  given  y  and  q) ,  with  test  values  clustered  in  particular 

around  zero  (parabolic  orbit)  and  close  to  64  (circular  orbit);  test 

values  for  x  also  clustered  about  zero  (perifocus),  with  a  maximum 
20 

value  of  10 

The  essentially  in-plane  tests,  just  described,  were  carried  out 
with  the  fixed  value  of  fcn  for  i  .  Tests  for  a  wide  range  of  i 
were  obviously  important,  however,  in  particular  with  values  clustered 
close  to  the  limiting  values  of  0  and  ir  ;  these  tests  were  carried 
out  with  a  varied  over  its  full  range,  but  x  (as  well  as  the  other 
quantities)  held  fixed.  Finally,  it  seemed  enough  to  carry  out  a 
limited  (non-systematic)  set  of  tests  with  the  values  of  ft  and  u> 
varied. 

As  already  indicated,  all  the  tests  met  the  goal  set  for  position 

error,  as  measured  by  |  Ar  |  /r  .  The  goal  was  also  met  for  velocity 

error,  as  measured  by  |  Af  |  /V  ,  in  most  of  the  tests,  but  it  was  not 

met  in  a  number  of  cases  involving  proximity  to  the  apofocus  of  a  near- 

rectilinear  orbit.  Since  V  is  very  small  in  such  circumstances,  the 

requirement  becomes  severe,  but  failure  to  achieve  it  can  be  ascribed  to 

a  specific  technical  point:  when  E  is  close  to  x  ,  it  is  impossible 

to  restrict  the  relative  error  in  sin  E  ;  but  V  is  proportional  to 

K 

sin  E  ,  so  when  r  is  dominated  by  (because  *  0)  the  loss  of 
(relative)  accuracy  in  r  is  inevitable;  another  way  of  putting  this 
point  is  that  5V/5x  is  large  in  comparison  with  V/x  ,  so  the  relative 
error  in  computing  x  (within  PV2ELS)  gets  magnified  in  the  subsequent 
return  to  velocity.  If  we  measure  relative  error  by  j  Air  j  /W  ,  how- 
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ever,  where  W  is  max(V,  /a),  not  just  V  ,  then  the  goal  of  2  x  10 


is  met  in  all  cases.  A  similar  problem  might  have  been  expected  with 
position  error  when  r  is  very  small;  but  for  values  of  E  close  to 
zero,  the  relative  error  in  sin  E  is  kept  low,  whilst  values  of  E 
close  to  2ir,  4t r  etc  cannot  arise  with  the  composite  function 
fff  ^(x)1  ,  since  f  ^  selects  for  minimum  t  as  explained  in 
Section  4. 


6.  On  the  Singularities 

Not  counting  the  rectilinear-orbit  singularity  (defined  by  q  ■  0  ,  but 
of  little  practical  interest),  there  are  three  singularities  that  have 
long  been  recognized  as  possible  sources  of  difficulty  in  computational 
analysis  based  on  the  usual  elliptic  elements:  the  singularity  at 
e  ■  0  ;  and  the  pair  of  singularities  at  i  •  0  and  i  ■  ir  .  It  is 
easy  to  transform  the  elements  so  as  to  eliminate  any  one  of  these 
singularities,  and  it  is  not  hard  to  find  a  transformation  that  sim¬ 
ultaneously  eliminates  the  eccentricity  singularity  and  one  of  the 
inclination  singularities  -  this  is  achieved  with  the  'equinoctial 
elements'  of  Broucke  and  Cefola  (1972),  for  example.  The  elimination  of 
all  three  singularities  together  is  much  more  difficult,  however, 
though  it  was  achieved  by  Cohen  and  Hubbard  (1962)  with  the  element  set 
defined  by 

qQ  -  p^  cos  ^i  cos  ^(£5  +  w  +  o)  ,  q^  *  p^  sin  ^1  cos  *s(£2  -  u  -  o)  , 

v  V 

q£  “  p  sin  4ji  sin  ^(£2  -  u  -  o)  ,  q^  ”  P  cos  ^i  sin  ^(£2  +  id  +  o)  , 

e  ■  e  cos  a  and  e  m  -  e  sin  o  , 

x  y  ’ 

where  o  is  the  mean  anomaly  at  epoch. 

The  trouble  with  the  partial  or  total  elimination  of  singularity, 
in  element  sets  such  as  the  ones  just  referred  to,  is  that  it  only 


applies  to  elliptic  orbits  and  is  incompatible  with  universal  appli¬ 
cability.  By  the  same  token,  there  is  no  simple  transformation  of  the 
universal  elements  recommended  here  that  eliminates  singularity  whilst 
retaining  universality.  We  have  already  seen,  on  the  other  hand,  that 
singularity  is  of  little  consequence  so  far  as  conversion  procedures  are 
concerned.  But  up  to  now  we  have  only  been  concerned  with  fixed  values 
of  the  relevant  elements;  when  variations  in  these  elements  are  con¬ 
sidered,  we  cannot  dismiss  the  singularity  question  so  lightly. 

Suppose  we  have  a  large  computer  program  involving  differential 
changes  in  the  values  of  a  set  of  orbital  elements  subject  to  singu¬ 
larity;  these  changes  may  arise,  in  particular,  from  the  evaluation  of 
perturbation  formulae,  or  from  the  use  of  observations  in  orbit  deter¬ 
mination.  If  the  program  is  used  with  an  element  set  in  the  vicinity  of 
a  singularity,  it  is  likely  that  unacceptable  error  will  result,  and  the 
first  time  this  occurs  It  may  seem  that  we  have  to  rewrite  the  program 
in  terms  of  a  different  (singularity  free)  element  set.  This  could  be  a 
daunting  prospect,  and  an  especially  unattractive  one  if  it  could  be 
fftreseen  that  the  efficiency  of  the  new  program  would  be  significantly 
reduced  by  the  complexity  of  the  new  elements.  However,  such  a  radical 
rewrite  is  almost  certainly  unnecessary,  as  it  should  be  possible  to 
maintain  the  original  elements  with  only  a  short-term  switch  to  non¬ 
singular  equivalents. 

To  illustrate  the  principle,  consider  just  the  singularity  associ¬ 
ated  with  e  ■  0  .  For  simplicity,  we  suppose  that  the  computer  program 
has  been  written  in  terms  of  the  usual  elliptic  elements,  rather  than 
our  universally  applicable  set,  but  this  is  a  point  of  minor  consider¬ 
ation.  Unless  e  is  exactly  zero,  we  have  been  able  to  compute  changes 
in  all  the  elements,  including  superficially  meaningless  values  of 
and  6M  .  Applying  these  changes  directly  could  lead  to  serious  error 
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in  e,  a)  and  M  ,  attributable  to  truncation  effects  in  the  underlying 
Taylor  expansions,  but  such  error  could  be  avoided  by  just  a  temporary 
switch  from  the  elements  e,  ui  and  M  to  n  and  U  ,  where 

(5.  h,  U)  ■  (e  cos  u),  e  sin  ui,  M  +  w); 

the  switch  is  implemented  on  the  basis  that  e'  ,  w '  and  M'  ,  the  new 
elements  we  require,  are  recoverable  from 


e '  cos  o> ' 

-  ?' 

-  Z  +  <5?  - 

(e  +  6e) 

cos  u  -  (e  6o>)  sin  u> 

e '  sin  in  ’ 

-  n' 

*  n  +  6n  * 

(e  +  5e) 

sin  u  +  (e  6o>)  cos  u> 

M'  +  «' 

-  O' 

-  U  +  6U  - 

(M  +  6M) 

+  (m  +  6w )  . 

Truncation  error  can  be  entirely  avoided  by  the  artifice  just 
illustrated  (except  perhaps  at  the  points  of  exact  singularity),  but 
rounding  error  can  still  be  a  problem.  Depending  on  the  application,  it 
should  be  possible  to  deal  with  this  via  minimal  software  modifications, 
in  particular  in  the  computation  of  partial  derivatives.  In  the  pertur¬ 
bation  context,  the  procedure  for  singularity  avoidance  has  been 
indicated  before,  by  Gooding  (1983). 

7.  Conclusion 

The  only  truly  'universal'  elements  are  perhaps  the  epochal  components 

of  position  and  velocity,  and  papers  such  as  Shepperd's  (1985)*  give 

*  Two  points  in  this  generally  excellent  paper  illustrate  the  weaknesses 
of  computing  procedures  that  are  almost  totally  blind  to  orbital  type. 
First,  Shepperd  recommends  the  use  of  a  certain  universal  variable,  u  , 
in  terms  of  which  an  'intermediate  parameter'  q  ,  where 
q  ■  ctu2/(l  +  au^)  ,  is  the  argument  of  a  particular  continued-fraction 
expansion,  Gs(q)  •  For  elliptic  orbits,  q  is  actually  sin2^AE  , 
where  AE  -  E  -  EQ  ,  so  that  qmax  *  °*5  »  as  Shepperd  remarks.  But  for 
this  value  of  q  ,  18  iterations  of  the  recursive  evaluation  of  G5<q) 
are  needed  to  give  14-decimal-digit  accuracy,  whereas  the  same  accuracy 
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elegant  algorithms  for  the  propagation  (or  transition)  of  these  from  one 
instant  of  time  to  another.  Only  one  element  of  the  traditional  type 
figures  in  such  algorithms,  namely,  the  energy-equivalent  inverse  of  the 
semi -major  axis.  However,  total  universality  is  not  available,  even 
with  an  approach  as  general  as  this,  since  the  generalized  anomaly  vari¬ 
able  has  to  be  range-restricted  when  an  elliptic  orbit  is  identified  and 
more  than  a  single  revolution  within  this  ellipse  is  involved.  Though 
text-books  increasingly  reflect  the  attractions  of  the  very  general 
approach,  they  also  continue  to  recognize  the  utility  of  traditional 
element  sets,  and  in  particular  of  sets  for  which  five  of  the  six 
elements  are  Independent  of  orbital  position.  It  has  been  the  main 
objective  of  the  present  paper  to  establish  that  the  most  familiar  of 
all  element  sets  is,  with  only  slight  modification,  of  universal 
application,  so  long  as  the  necessary  conversion  procedures  (to  and  from 
position  and  velocity)  are  carefully  programmed.  In  particular,  optimum 
numerical  accuracy  can  only  be  maintained  if  the  procedures  respect  the 
necessity  for  different  types  of  formulae  to  be  employed  Internally, 
according  to  the  type  of  orbit.  , 

The  universality  objective  (for  an  element  set)  Is  not  compatible 
with  freedom  from  singularity,  but  in  the  majority  of  applications 
singularity  is  either  of  little  consequence  or  can  be  dealt  with  easily. 
Hence  the  existence  of  singularities  for  a  universally  applicable 
element  set  should  not  be  regarded  as  a  major  defect. 


can  be  obtained  with  much  greater  efficiency  if  the  G5  function  Is 
evaluated,  in  a  specifically  elliptic  formulation,  as 
15(6AE  -  8  sin  AE  +  sin  2AE)/96  sin^AE  .  The  other  point  is  more 
serious,  in  that  Shepperd  suggests  that  use  of  u  as  argument  of  a 
universal  Kepler’s  equation  automatically  eliminates  the 
s low-convergence  problem  noted,  in  particular,  by  Odell  and  Gooding 
(1986).  This  is  a  misconception,  however;  with  Newton-Raphson 
iteration,  the  problem  can  only  be  universally  eliminated  by  careful 
choice  of  a  starting  formula. 
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Four  Fortran-77  subroutines,  which  implement  the  conversion  procedures 
described  in  the  paper,  are  listed  in  Appendices  A  to  D. 


oooooooo 
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Appendix  A 

SUBROUTINE  ELS2PV 

SUBROUTINE  ELS2PV  (GM,  AL,  Q,  OM,  TAU,  R,  U,  VR,  VT) 

ALGORITHM  FOR  TWO-DIMENSIONAL  CONVERSION 
FROM  ORBITAL  ELEMENTS  TO  POSITION  AND  VELOCITY. 

INPUT  ARGUMENTS  ARE:  GM  (G*M) ,  AL(PHA)  (GM/A) , 

Q  (PERI  DISTANCE),  OM (EGA)  (ARG-PERI  RELATIVE  TO 
ASSUMED  REFERENCE  DIRECTION)  AND  TAU  (TIME  FROM  PERI). 
OUTPUT  ARGUMENTS  ARE:  R  (RADIAL  DISTANCE) , 

U  (ANGLE  FROM  REFERENCE  DIRECTION) ,  VR  (RADIAL  VELOCIY) 
AND  VT  (TRANSVERSE  VELOCITY:  .GE.O). 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

PARAMETER  (PI  *  3 . 14159265358979323846264338328D0 , 

1  FOURPI  =  4D0*PI) 

IF  (AL.EQ.ODO)  THEN 

C  (PARABOLA  -  GM  CANNOT  BE  ZERO) 

D  «  DCBSOL(0 . 5D0/GM,  Q,  1. 5D0*GM*TAU) 

R  *  Q  +  0 . 5D0*D*D/GM 
H  =  DSQRT (2D0*GM*Q) 

V  -  2D0*DATAN2 (D,  H) 

ELSE 

C  (ELLIPSE  OR  HYPERBOLA) 

El  =  AL*Q 
E  -  GM  -  El 
EP1  »  GM  +  E 
H  -  DSQRT (Q*EP1) 

ALP  *  DABS (AL) 

RTAL  =  DSQRT (ALP) 

C  (LAST  6  ITEMS  COULD  BE  SAVED  IF  REPEATING  GM,  AL  &  Q) 

EM  =  TAU  * ALP  *RTAL 
IF  (AL.GT.ODO)  THEN 

C  (ELLIPSE  -  GM  CANNOT  BE  ZERO) 

EE2  -  0 . 5D0*EKEPL(EM/GM,  El/GM) 

S2  =  DSIN (EE2 ) 

C2  »  DCOS (EE2) 

R  =  Q  +  2D0*E*S2*S2/AL 
D  -  2D0*E*S2*C2/RTAL 

V  =  2D0*DATAN2 (EP1*S2 ,  H*RTAL*C2 ) 

EMV  -  EM/GM  -  V 

V  ■  V  +  FOURPI*DSIGN (DINT (DABS (EMV/FOURPI)  +  0.5D0),  EMV) 
ELSE 

C  (HYPERBOLA) 

S  =  SHKEPL(EM/E,  -El/E) 

S2  =  S*S 

C  =  DSQRT (IDO  +  S2) 

S2  =  S2/(C  +  IDO) 

R  =  Q  -  E*S2/AL 
D  =  E*S/RTAL 

V  *  DATAN2 (S*H*RTAL,  -GM*S2  -  El) 

END  IF 

END  IF 

C  (ALL  ORBITS) 

U  =  OM  +  V 
VR  =>  D/R 
VT  =  H/R 
RETURN 
END 


o  o 
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Appendix  B 
SUBROUTINE  ELS3PV 

SUBROUTINE  ELS3PV  (GM,  AL,  Q,  El,  BOM,  OM,  TAU, 

1  X,  Y ,  Z,  XDOT,  YDOT,  ZDOT) 

ALGORITHM  FOR  THREE-DIMENSIONAL  CONVERSION 
FROM  ORBITAL  ELEMENTS  TO  POSITION  AND  VELOCITY. 
IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

CALL  ELS2PV  (GM,  AL,  Q,  OM,  TAU,  R,  U,  VR,  VT) 

C  =  DCOS(U) 

S  =  DSIN(U) 

XI  *  R*C 
Y1  =  R*S 

X2  =  VR*C  -  VT*S 
Y2  =  VR*S  +  VT*C 
C  =  DCOS(EI) 

S  =  DSIN(EI) 

2  =  Y1*S 
VI  *  Y1*C 
ZDOT  -  Y2*S 
Y2  *  Y2*C 
C  -  DCOS(BOM) 

S  -  DSIN(BOM) 

X  -  X1*C  -  Y1*S 
Y  *  X1*S  +  Y1*C 
XDOT  *  X2*C  -  Y2*S 
YDOT  =  X2*S  +  Y2*C 
RETURN 
END 
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Appendix  C 
SUBROUTINE  PV2ELS 

SUBROUTINE  PV2ELS  (GM,  R,  U,  VR,  VT,  ALf  Q,  OM,  TAU) 

C  ALGORITHM  FOR  TWO-DIMENSIONAL  CONVERSION 

C  FROM  POSITION  AND  VELOCITY  TO  ORBITAL  ELEMENTS. 

C  INPUT  ARGUMENTS  ARE:  GM  (G*M) ,  R  (RADIAL  DISTANCE) , 

C  U  (ANGLE  FROM  ASSUMED  REFERENCE  DIRECTION) , 

C  VR  (RADIAL  VELOCITY)  AND  VT  (TRANSVERSE  VELOCITY:  .GE.O). 

C  OUTPUT  ARGUMENTS  ARE:  AL(PHA)  (GM/A) ,  Q  (PERI  DISTANCE), 

C  OM(EGA)  ( ARG-PERI  RELATIVE  TO  REFERENCE  DIRECTION) 

C  AND  TAU  (TIME  FROM  PERI). 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

LOGICAL  L 

PARAMETER  (PI  *  3 . 14159265358979323846264338328D0 , 

1  TWOPI  -  2D0*PI ,  SW  =  0.25D0,  L  =  .FALSE.) 

C  (ALL  ORBITS) 

VSQ  *  VR*VR  +  VT*VT 
AL  *  2D0*GM/R  -  VSQ 
ALP  *  DABS (AL) 

RTAL  -  DSQRT (ALP) 

D  =  R*VR 

H  «  R*VT 

P  =  H*H 

ESQ1  -  P*AL 

ES  »  D*RTAL 

ESES  *  ES*ES 

EC  -  R*VSQ  -  GM 

ECEC  *  EC*EC 

IF  (AL.GT.ODO)  THEN 

C  (ONE  ESQ  FORMULA  SUPERIOR  FOR  THE  ELLIPSE) 

ESQ  »  ECEC  +  ESES 
ELSE 

C  (DIFFERENT  FORMULA  SUPERIOR  FOR  THE  HYPERBOLA) 

ESQ  «  GM*GM  -  ESQ1 
END  IF 

E  »  DSQRT (ESQ) 

Q  *  P/  (GM  +  E) 

IF  (AL. EQ. 0D0)  THEN 
C  ( PARABOLA) 

TAU  *  D*(2D0*Q  +  R)/(3D0*GM) 

V  »  2DO*DATAN2(VR,  VT) 

ELSE  IF  (E.EQ.ODO)  THEN 

C  (CIRCLE) 

TAU  *  ODO 

V  *  ODO 
Ea"E 
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(continuation  of  PV2ELS ) 

C  (ELLIPSE  OR  HYPERBOLA) 

El  =  AL*Q 

IF  (AL.GT.ODO)  THEN 
C  (ELLIPSE) 

EH  =  DAT  AN  2  (ES,  EC) 

IF  (GM*EH*EH/6D0  +  El  .GE.  GM*SW)  THEN 
C  (GENERAL  CASE) 

EM  =  GM*£H  -  ES 
ECESQ  =  GM*EC  -  ESQ 
ELSE 

C  (FOR  El  &  EH  BOTH  NEAR  ZERO) 

EM  =  GM*EMKEP (El/GM,  EH) 

ECESQ  =  (ESQ1*ECEC  -  ESQ*ESES)/ (ESQ  +  GM*EC) 

END  IF 
ELSE 

C  (HYPERBOLA) 

EH  =  DASINH (ES/E) 

IF  (GM*EH*EH/6D0  -  El  .GE.  GM*SW)  THEN 
C  (GENERAL  CASE) 

EM  =  ES  -  GM*EH 
ECESQ  **  ESQ  -  GM*EC 
ELSE 

C  (FOR  El  &  EH  BOTH  NEAR  ZERO) 

EM  =  E*SHMKEP (-E1/E,  ES/E) 

ECESQ  *  - (ESQ1*ECEC  +  ESQ*ESES) / (ESQ  +  GM*EC) 

END  IF 
END  IF 

C  (ELLIPSE  OR  HYPERBOLA  STILL) 

EN  =  ALP*RTAL 
TAU  *  EM/EN 

V  -  DATAN2(ES*H*RTAL,  ECESQ) 

END  IF 

C  (ALL  ORBITS) 

OM  =  U  -  V 

IF  (L  .AND.  AL.GT.ODO)  THEN 

C  (FOR  ELLIPSE,  ADJUST  REVOLUTIONS  IF  REQUIRED  (USING  L) ) 

ADJ  =  TWOPI *DSIGN( DINT ( DABS (OM/TWOPI)  +  0.5D0),  OM) 

OM  =  OM  -  ADJ 
TAU  *  TAU  +  ADJ/EN 
END  IF 
RETURN 
END 
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Appendix  D 
SUBROUTINE  PV3ELS 

SUBROUTINE  PV3ELS  (GM,  X,  Y,  Z,  XDOT,  YDOT,  ZEX3T, 

1  AL,  Q,  El,  BOM,  OM,  TAU) 

C  ALGORITHM  FOR  THREE-DIMENSIONAL  CONVERSION 

C  FROM  POSITION  AND  VELOCITY  TO  ORBITAL  ELEMENTS. 

IMPLICIT  DOUBLE  PRECISION  (A-H,0-Z) 

PARAMETER  (PI=3 . 14159265358979323846264338328D0 ,  HALFPI=PI/2D0) 
XSQYSQ  =  X*X  +  Y*Y 
RSQ  =  XSQYSQ  +  Z*Z 
R  =  DSQRT(RSQ) 

VR  =  (X*XDOT  +  Y*YDOT  +  Z*ZDOT)/R 
HX  =  Y*ZDOT  -  Z*YDOT 
HY  =  Z*XDOT  -  X*ZDOT 
HZ  =  X*YDOT  -  Y*XDOT 
HSQ  =  HX*HX  +  HY*HY  +  HZ*HZ 
IF  (HSQ.EQ.ODO)  THEN 
C  (RECTILINEAR  ORBIT) 

El  *  HALFPI 

IF  (XSQYSQ. EQ .0D0)  THEN 

C  (AXIAL  ORBIT) 

BOM  =  0D0 
ELSE 

C  (GENERAL  RECTILINEAR  ORBIT) 

BOM  =  DATAN2 (Y,  X) 

END  IF 

U  =  DATAN2 (Z ,  DSQRT (XSQYSQ) ) 

VT  =  0D0 
ELSE 

C  (NON-DEGENERATE  ORBIT) 

BX  =  HY*Z  -  HZ*Y 
BY  =  HZ*X  -  HX*Z 
BZ  *  HX*Y  -  HY*X 
HX  =  Y*BZ  -  Z*BY 
HY  =  Z*BX  -  X*BZ 
HZ  =  X*BY  -  Y*BX 
W  =  HX*HX  +  HY*HY 
H  =  DSQRT (W  +  HZ*HZ) 

El  =  DATAN2 (DSQRT (W) ,  HZ) 

IF  (W.EQ.ODO)  THEN 

C  (ORBIT  IN  REFERENCE  PLANE) 

BOM  =  0D0 

U  =  DATAN2 (Y*DSIGN ( IDO ,HZ) ,  X) 

ELSE 

C  (GENERAL  ORBIT) 

BOM  =  DATAN2 (HX,  -HY) 

U  =  DATAN2 (H*Z ,  RSQ*BZ) 

END  IF 

VT  =  H/(R*RSQ) 

END  IF 

CALL  PV2ELS  (GM,  R,  U,  VR,  VT,  AL,  Q,  OM,  TAU) 

RETURN 

END 
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